NASA'. NASA-TP-2385 19850007385 

Technical 

Paper 

2385 

January 1 985 # 

A Comparative Study 
of the Nonuniqueness 
Problem of the 
Potential Equation 

M. D. Salas, 

A. Jameson, 
and R. E. Melnik 



SEP . 1 » 

LANGLEY RESEARCH CENTER 
i JBRARY NASA, HAMPTON, VA. 




3 1 176 01331 7756 




NASA 

Technical 

Paper 

2385 

1985 


A Comparative Study 
of the Nonuniqueness 
Problem of the 
Potential Equation 


M. D. Salas 

Langley Research Center 
Hampton, Virginia 

A. Jameson 

Princeton University 
Princeton, New Jersey 

R. E. Melnik 

Grumman Aerospace Corporation 
Bethpage, New York 


NASA 

National Aeronautics 
and Space Administration 


Scientific and Technical 
Information Branch 



Summary 

The nonuniqueness problem occurring at transonic 
speeds with the conservative potential equation is in- 
vestigated numerically. Evidence is given supporting 
the thesis that the nonuniqueness problem is inherent 
in the potential differential equation rather than a con- 
sequence of the finite-difference approximation. Results 
are presented from an extensive comparative study be- 
tween potential and Euler calculations for flow past 
two-dimensional airfoil profiles. This study indicates 
that the nonuniqueness problem is not an inviscid phe- 
nomenon, but results from the approximate treatment 
of shock waves inherent in the conservative potential 
model. A more restrictive bound on the limit of valid- 
ity of the conservative potential model is suggested. 

Introduction 

When the tangential stresses acting on a volume of 
fluid are small compared with the normal pressures ex- 
erted by the fluid, the fluid may be assumed to be invis- 
cid. The mathematical expressions describing the con- 
servation of mass, momentum, and energy of an inviscid 
flow are known as the Euler equations. In the absence 
of rotational effects, the inviscid flow can be described 
by a single conservation law, the conservation of mass, 
which can be expressed in terms of a function <3>, the 
velocity potential. The irrotational assumption is not 
always valid. Rotational effects can occur in an inviscid 
flow as a result of variations in free-stream stagnation 
conditions, because of the presence of shock waves, or 
because of satisfying boundary conditions, as for ex- 
ample in the case of a lifting delta wing with a sharp 
leading edge. For flow past two-dimensional airfoils im- 
mersed in a uniform upstream far field, rotational ef- 
fects are introduced at shock waves. Strictly speaking, 
within the inviscid assumption, these flows are governed 
by the Euler equations. However, if shock waves are suf- 
ficiently weak, the irrotational assumption should pro- 
vide a reasonable approximation of the flow field. Shock 
waves, under this approximation, are replaced by sharp 
isentropic recompressions which approximate the true 
Rankine-Hugoniot shock up to second order in terms of 
the shock wave strength. 

Under the irrotational assumption, the conservation 
of mass equation can be written in quasi-linear or 
nonconservative form as 

(a 2 - $ 2 )$ I:I - + (a 2 - <S> 2 y )<b yy = 0 

a = a(<D 2 + d> 2 ) 

or in conservation form as 

(f&x)x + {P$y)y - 0 1 

p = p(*l + Ql) } 


where a is the speed of sound, p is the density, $ is 
the potential, and the subscripts denote partial deriva- 
tives. In the absence of shock waves, both forms of the 
equation are equivalent. When shock waves are present, 
use of the conservative form in conjunction with an “en- 
tropy” condition guarantees (see ref. 1) that the proper 1 
shock jump condition is satisfied. On the other hand, 
use of the nonconservative form leads to an ambigu- 
ous jump across the shock waves. With few exceptions, 
present-day state-of-the-art working transonic potential 
equation codes use the conservative form. 

For subcritical flows, Bers (ref. 2) has shown that 
the solution for flow past an airfoil section satisfying 
the Kutta condition is unique. In the transonic regime, 
little has been proven rigorously in terms of uniqueness. 
However, in a recent paper (ref. 3), Steinhoff and Jame- 
son showed that numerical solutions of the full potential 
equation in conservative form for flow past an airfoil are 
not unique. In the same paper, a meticulous study was 
conducted which indicated that the anomaly is inherent 
in the partial differential equation and not a result of 
its discrete representation. Having reached this conclu- 
sion, the authors went on to conjecture the possibility 
of the anomaly being (1) an inviscid phenomenon, per- 
haps inherent to the more exact Euler equations; and 
(2) perhaps a contributing factor in the occurrence of 
buffeting. 

Given the preeminent position today of the poten- 
tial approximation in the transonic range for the design 
and analysis of airfoil profiles, as documented in refer- 
ence 4, it is unnecessary to dwell on the importance of 
investigating the conjectures made in reference 3. Such 
an investigation is the purpose of this paper. The paper 
is divided into three main parts. First, a brief review of 
Steinhoff and Jameson’s study and some new findings 
made in the course of this work are presented. Second, 
the Euler code described in reference 5 is introduced, 
and both improvements and the validation procedure 
used to determine its accuracy are described. Finally, 
results from a comparative study of the nonuniqueness 
problem for the potential and Euler models are pre- 
sented and discussed. 

This paper provides results which are more compre- 
hensive than the preliminary results reported in refer- 
ence 6 and corrects some errors in an earlier version by 
the same title, reference 7. 

Symbols 

a local speed of sound 

C p pressure coefficient 


1 By proper, we mean consistent with the mathematical model 
being solved. In this case, the model requires the conservation of 


( 2 ) 


mass. 



c airfoil chord 

Cd section drag coefficient 

ci section lift coefficient 

c m section pitching-moment coefficient 

D damping operator defined by equa- 

tion (6) 

Dx damping operator defined by equa- 

tion (7) 

Dy damping operator defined by equa- 

tion (8) 

d term defined by equation (9) 

h cell area 

Jmax number of cells in Y-direction 

L wind-tunnel half-height 

M Mach number 

r radial distance measured from center of 

airfoil 

S similarity entropy 

t time 

u component of velocity in x-direction 

v component of velocity in y-direction 

w vector of conservative unknowns 

X computational coordinate along ring of 

“O” mesh 

x Cartesian coordinate 

Y computational coordinate normal to ring 

of “O” mesh 

y Cartesian coordinate 

a angle of attack 

P compressibility factor, (1 — M^) 1 / 2 

T circulation 

7 ratio of specific heats 

A increment 

e second-order artificial damping 
coefficient 

e fourth-order artificial damping 
coefficient 

9 polar angle measured counterclockwise 

from trailing edge 

k term defined by equation (13) 


p density 

r ratio of airfoil thickness to airfoil chord 

$ velocity potential 

(f) reduced velocity potential, 

0 = $ — ( a 00 M 00 )(xcosa + y sin a) 

X transonic similarity, 

(l-M^)/[r( 7 +l)M^] 2 /3 

Subscripts: 

d value downstream from shock wave 

F far-field value 

i cell index in X-direction 

j cell index in Y-direction 

o value at which nonuniqueness occurs 

s value upstream from shock wave 

oo free-stream value 

Review of the Potential Anomaly 

While introducing lifting capability to the full- 
potential, conservative, nonlifting, multigrid code de- 
scribed by Jameson in reference 8, Steinhoff and Jame- 
son (ref. 3) discovered that three converged solutions 
could be obtained for a symmetric airfoil at zero angle 
of attack and a fixed free-stream Mach number. These 
solutions consisted of a zero-lift symmetric solution and 
two mirror- image asymmetric solutions with large ab- 
solute levels of lift. (See figs. 1 and 2 of ref. 3.) The 
latter two solutions could be obtained by perturbing 
the zero-lift solution, which was unstable. The multi- 
grid capability of the code made it feasible to converge 
all three solutions to machine accuracy. 2 This unques- 
tionably demonstrated the multiplicity of solutions for 
the discrete system of equations. To determine if the 
anomaly was characteristic of the partial differential 
equation being modeled, the following tests were made 
in reference 3 with FLO 36: 

1. The solutions were checked to determine if they 
had the expected far-field decay corresponding to a 
vortex and a doublet in a uniform stream (ref. 9). They 
did. In addition, the effect of the location of the outer 
boundary was investigated. It affected the lift level 
slightly, but the nonuniqueness problem persisted. 

2. Originally, the anomaly was observed for 
an airfoil with nonzero trailing-edge included angle. 
To determine the effect of imposing the Kutta con- 
dition, a Joukowski profile was investigated. The 

2 Residuals were of the order of 10“ 12 in a CDC® CYBER 
203 computer for the cases reported here. 


2 



Joukowski profile, with its cusped trailing edge, requires 
a smooth matching of the pressures on the upper and 
lower surfaces. The numerical solutions showed the 
proper trailing-edge behavior, but again the anomaly 
persisted. 

3. The convergence of these solutions with mesh 
refinement was investigated by using grids with 96 x 24, 
128 x 32, 192 x 48, 256 x 64, and 384 x 192 cells in the 
X- and Y-directions, respectively. All solutions were 
found to change little with mesh refinement. 

4. The effect of artificial viscosity was tested by 
using both first- and second-order accurate expressions 
for this term, with no effect on the anomaly. 

5. Since all calculations were originally made with 
an u O” mesh defined by a conformal circle mapping, 
calculations were tried with a “C” mesh. The multiple 
solutions were also observed with the “C” mesh. 

6. Similar results were obtained with the conserva- 
tive finite-difference, potential code described in refer- 
ence 10. 

In addition, during the course of the present inves- 
tigation, we found the following results: 

1. The multiple solutions also occurred with the 
conservative potential code described in reference 11 
(TAIR). 

2. A standard feature of these potential codes is 
the requirement that the reduced far-field potential 

satisfy the leading term of the compressible vortex 
solution; namely, 

r 

4>p = — tan -1 [ft tan(0 — a)] (3) 

0=(1 -Ml) 1 ' 2 (4) 

where 6 is the polar angle, and T is the value of the 
circulation that satisfied the Kutta condition. Because 
T evolved as part of the solution, it is reasonable to 
be suspicious of this boundary condition as a possible 
cause of the anomaly. To check this possibility, calcula- 
tions were made with the “wind-tunnel” code described 
in reference 12. In this code, the far-field boundary con- 
dition corresponds to no flow through the tunnel walls; 
that is, 


Essentially “free- air” solutions are generated with L — 
50c. The multiple solutions were also found with this 
code. 

3. The same anomalous behavior was found with 
the conservative, Cartesian coordinate code described 
in reference 13. 

4. No multiple solutions were found with the non- 
conservative code described in reference 14 (FLO 12). 

Although the above tests are not conclusive, they 
make a strong case for the argument that the problem 


exists at the differential equation level of the conserva- 
tive potential formulation. Therefore, with no evidence 
to the contrary, we will proceed under that assumption. 

An interesting feature of the results presented in ref- 
erence 3 was the occurrence of a “gap” in the lift curve. 
(See fig. 1.) This gap was explained in reference 3 in 
terms of a hysteresis effect. In the present investigation, 
it was determined that the gap comes about because the 
lift becomes a multivalued function of the angle of at- 
tack in this range. However, since the angle of attack 
remains a single- valued function of the lift, the prob- 
lem can be made well posed by prescribing the lift and 
letting the numerical solution determine the angle of 
attack which satisfies the Kutta condition. (A similar 
technique was successfully applied by the principal au- 
thor to another nonunique problem, ref. 15.) In this 
manner, the complete lift curve for an NACA 0012 air- 
foil has been evaluated at = 0.83 with FLO 36, a 
modified version of the code described in reference 3. 
The results are shown in figure 2. 

For the NACA 0012 airfoil, the results presented 
in figure 2 show the anomaly occurring in the neigh- 
borhood of the zero angle of attack. It should not be 
inferred from this that the anomaly occurs only at this 
Mach number and narrow range of angles of attack. 
At any given Mach number in the transonic regime, 
a range of angles of attack can be found where the 
nonuniqueness occurs. For example, for the NACA 0012 
airfoil at M 0 0 — 0.79, it appears at around a = 1°; at 
Moo = 0.60, it appears at around a = 9°. 

It has been speculated that the anomaly is a result 
of an interaction between an upper-surface shock and 
a lower-surface shock. However, the anomaly has been 
observed in a number of cases in which only a single 
shock (either on the lower or upper surface) is present. 
A case in point is the flow past an RAE 2822 profile 
at Mqo = 0.75 and a = 1°. The lift curve for these 
conditions is shown in figure 3. The single shock wave 
occurs on the upper surface, as shown in figure 4. 

Nonuniqueness is not new to potential theory. An 
incompressible flow past a lifting airfoil with a sharp 
trailing edge has a one-parameter family of solutions. In 
this case, the relevant physical solution is singled out by 
the Kutta condition. At supercritical speeds, solutions 
with “expansion” shocks are mathematically valid solu- 
tions of the governing differential equation. Expansion 
shocks are ruled out in favor of the physically relevant 
compression shock by invoking the Second Law of Ther- 
modynamics. In fact, these constraints are now consid- 
ered integral parts of a numerical algorithm. Why not 
resolve the present dilemma by imposing some new con- 
straint? If we consider figure 2 again, at zero angle of 
attack, the symmetric solution corresponding to zero 
lift seems to be the solution that we would want to sin- 
gle out as relevant with the new constraint. However, 
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consider what happens if the angle of attack is increased 
by a small amount. Which of these solutions is physi- 
cally relevant? The two with the largest absolute levels 
of lift can be ruled out on the grounds that they require 
a discontinuous behavior as a goes to 0°. However, the 
remaining solution does not seem physically relevant ei- 
ther, since it predicts a negative lift-curve slope. The 
nonuniqueness problem we are facing is unique in that 
none of the solutions available seem to be physically rel- 
evant! Indeed, if a new constraint is found, its role will 
not be to single out one of the three solutions presently 
available, but to find a new solution. 

The question now is, Is this anomaly a problem of 
the conservative potential approximation in the tran- 
sonic range or a problem of the inviscid flow? To in- 
vestigate this question, we propose to do a systematic 
search with an Euler flow code. However, before pro- 
ceeding to this step, it is necessary to establish the va- 
lidity of the numerical results of the Euler code. This 
is the next step in our investigation. 

Improvements and Validation of FLO 52-S 

The finite-volume, four-step Runge-Kutta code de- 
scribed in reference 5 (FLO 52) was chosen for this 
study because of its fast convergence rate. The code 
solves the time-dependent continuity, x-momentum, y- 
momentum, and energy equations in conservative form. 
To prevent typical central-difference, odd-even oscilla- 
tions and to capture shock waves, a blended fourth- and 
second-order explicit artificial damping term is added 
to the equations. To accelerate the rate of convergence, 
each mesh cell is advanced at its own local time step, 
A £, and the governing equations are modified by an en- 
thalpy damping term. The calculations were done with 
a nearly orthogonal “O” mesh with an inverse radial 
transformation to cluster mesh points near the airfoil. 
A detailed account of the discretization, damping, and 
acceleration technique is given in reference 5. 

In its original version, the code consistently under- 
predicted the lift coefficient from potential calculations 
by as much as 10 percent for lifting subcritical flows. 
In addition, isomach plots revealed a disturbing behav- 
ior near the airfoil surface. As the isomach lines ap- 
proached the surface, they turned abruptly. This turn- 
ing indicated an incorrect “boundary-layerlike” behav- 
ior. This was also evident in flow-field plots of entropy 
and total pressure losses. 

The cause of the problem was found to be associ- 
ated with the evaluation of the damping terms near the 
airfoil surface and in the far field. As described in ref- 
erence 5, the damping term is made up of contributions 
from the two coordinate directions X and Y with as- 
sociated indexes i and j, respectively. Let X be the 
coordinate along a ring of the “O” mesh and Y be the 


coordinate in the normal direction. The damping term 
is given by 

Dw = Dxw - h Dyw (6) 

D x w = d i+1/2tj - di-i /2 j (7) 

Dyw = d itj+1/2 - d itj _ i /2 (8) 

where w is the vector of conservative unknowns, and a 
typical d term is defined by 

A hid +1/2 L(2) a 

d iJ+l/2 = [eij+1/2 *V>ij+l/2 

~ £ i,j+l/2(^ W iJ+3/2 ~~ 2 AWij+x/2 

+ A^j_i/ 2 )] (9) 

The symbol h represents the cell area, e ^ and are 
coefficients, defined in reference 5, which depend on the 
second derivative of pressure, and 

AWij+1/2 — Wij+1 — Wi y j (10) 

The evaluation of the d terms presents some problems 
near the airfoil surface, j < 2, and near the last ring 
of the “O” mesh, j > / max - 1, because of the lack 
of information near these cells. In the original code, 
A w terms inside the airfoil and outside the last ring 
of the “O” mesh were constructed with a third-order 
extrapolation from inside the flow field. 

The “boundary-layerlike” behavior was eliminated 
by the following procedure. At a ghost point j = 
0 inside the airfoil, the pressure, density, and total 
enthalpy were obtained by linear extrapolation from the 
flow field. The velocity components at the cell center 
immediately next to the airfoil surface were decomposed 
into components normal and tangent to the surface. 
These were then reflected to obtain components at the 
ghost cell. With this information, it was possible to 
evaluate 

= w iy i - Wi^o 

The missing term A ^,_ 1 /2 was evaluated by equating 
it to Aw^i/ 2 . Although this procedure eliminated the 
“boundary-layerlike” behavior, it had little effect on 
the low-lift level being predicted. Surprisingly, the 
problem with the lift was corrected by modifying the 
evaluation of the far-field damping terms. The finite- 
volume integration is performed only up to J max — 1; 
values of w at J ma x are obtained by satisfying the 
far-field boundary conditions, which use no damping. 
Therefore, to evaluate the damping terms at J max — 1, 
only additional information for the term A^ ( j max+1 / 2 
is required. When this term was approximated by 

Aw i,J m * x + 1/2 = AWi,4„- 1/2 
the error in the lift level was reduced to about 1 percent. 
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The far-field boundary condition evaluation in the 
original code used the method of Rudy and Strikwerda 
(ref, 16), which introduces an additional free param- 
eter into the calculation. This method was replaced 
by the following procedure, which proved more robust 
and made a slight improvement in the convergence rate. 
At a subsonic point, a frame of reference perpendicu- 
lar and tangential to the last ring of the u O” mesh is 
constructed. In the plane made by the perpendicular 
to the ring of the “0” mesh and the time axis, the 
Riemann “invariant” that propagates along the char- 
acteristic coming into the computational region is pre- 
scribed from known free-stream values. The Riemann 
“invariant” propagated by the outgoing characteristic 
is extrapolated from computed values. From these two 
quantities, the speed of sound and the velocity com- 
ponent normal to the ring of the “O” mesh are deter- 
mined. If the point corresponds to an inflow point, the 
entropy and velocity components tangent to the ring of 
the “O” mesh are prescribed from free-stream values. If 
the point corresponds to an outflow point, the entropy 
and velocity components tangent to the ring of the “O” 
mesh are extrapolated from computed values. From 
these four quantities, the conservative unknown vector 
w is evaluated. For supersonic inflow, all quantities are 
prescribed, but for supersonic outflow, all quantities are 
extrapolated from computed values. 

The code was also modified to include, as an option, 
the effect of a far-field vortex. To this end, the circula- 
tion is approximated from the lift calculated by the in- 
tegration of the surface pressure. Then, the free-stream 
Cartesian velocity components modified to account for 
the vortex are defined as 


where 


u Q o — u Q o — k, sin 6 

(ii) 

Voo = V co + K, COS 0 

(12) 

27rr[l-M2 > sin 2 (0-c*)] 

(13) 

r = ^ ( -’&oo 

(14) 


and r, 0 are polar coordinates in the physical plane with 
the origin at the center of the airfoil. The free-stream 
thermodynamic variables are recomputed by using the 
magnitude of the modified free-stream velocity given 
by equations (11) and (12), the steady-state Bernoulli 
equation, and the known value of entropy at the point 
in question. In general, the far-field vortex had very 
little effect on the calculated results. This is probably 
the result of the location of the last ring of the “O” 
mesh, which in these calculations corresponded to about 
100 chords from the airfoil. 

The last major modification to the Euler code con- 


sisted of a second-order accurate integration of the nor- 
mal pressure gradient at the surface of the airfoil to eval- 
uate the surface pressure. As discussed in reference 5, 
because of the finite-volume formulation and the fact 
that the fluxes convected across the airfoil surface are 
zero, only the pressure needs to be evaluated at the sur- 
face. In this connection, it should be stressed that the 
evaluation of ghost points, discussed in relation to the 
damping terms, is not related to the evaluation of the 
surface boundary conditions. 

For reference purposes, we will designate the mod- 
ified Euler code as FLO 52-S. A large number of lift- 
ing subcritical cases were calculated with this code by 
using a 120 x 34 mesh and with the potential code 
FLO 36 by using a 192 x 32 mesh. Results are pre- 
sented for one of the “most difficult” subcritical lifting 
cases computed— an NACA 0012 airfoil at = 0.30 
and a. — 10°. A partial view of the mesh corresponding 
to 120 x 34 points is shown in figure 5. Figure 6 shows 
that second-order accuracy in the lift and drag coef- 
ficients is attainable with mesh refinement. Figure 7 
shows a typical convergence behavior. The residuals 
plotted in figure 7(a) are the root- mean-square devia- 
tions of Apj At, Apu/At, and Apvj At evaluated over 
the entire field. Figure 7(b) shows the convergence of 
lift, drag, and moment coefficients, each evaluated from 
the integration of the surface pressure. A comparison 
between FLO 36 and FLO 52-S for the surface pressure 
distribution, isomach contours, and streamline trajecto- 
ries is shown in figures 8 to 10. Practically no difference 
can be observed if the results of the two calculations are 
overlaid; we emphasize the results presented correspond 
to a case that showed the most disagreement with a po- 
tential calculation. 

At supercritical speeds, the potential calculations 
are only approximate solutions to the inviscid flow. 
Typically, these calculations have a stronger shock 
wave, which appears further aft than in Euler calcu- 
lations. Two typical results at = 0.82 and 0.86 
(both at zero angle of attack) are shown in figures 11 
and 12. A check on how well the shock wave is captured 
is shown in figure 13, where the calculated shock jump 
is compared with the exact jump. Taking into consid- 
eration the difficulty in reading the results because of 
the shock-wave smearing and the Zierep singularity, the 
agreement is quite good. For comparison, the results of 
the nonconservative potential code, FLO 12, are also 
included. 

Nonuniqueness Study 

The lift curve shown in figure 2 is repeated in fig- 
ure 14, along with the results from the Euler calcu- 
lation (FLO 52-S) and the nonconservative potential 
code FLO 12. As can be seen from this figure, the Eu- 
ler results do not show any anomaly at this particular 
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Mach number. Can it be that the Euler solution will 
be anomalous at some other Mach number? 

In order to study this, consider the angle made by 
the lift curve with the abscissa at zero lift. If this angle 
exceeds 90°, then we are into the anomalous region, 
since a positive angle of attack produces negative lift. 
Figure 15 shows the behavior of this angle as a function 
of free-stream Mach number for conservative potential, 
nonconservative potential, and Euler calculations. The 
figure shows that the conservative potential solution 
becomes progressively worse as the free-stream Mach 
number increases and eventually crosses over into the 
multiple solution region at about M ^ = 0.82. Both 
Euler and nonconservative potential solutions remain 
well-behaved throughout the Mach number range. The 
anomaly, therefore, is indicative of a breakdown of 
the conservative potential formulation, not of inviscid 
theory, per se. 

In figure 16, the plot of figure 15 is repeated in terms 
of transonic similarity parameters (ref. 17). In addition 
to the results obtained for an NACA 0012 airfoil, re- 
sults for an NACA 0006 airfoil are also shown in this 
figure. It is interesting to see that the conservative po- 
tential results appear to collapse into a single curve in 
the similarity plane. This is another indication that 
the numerical code is predicting the behavior expected 
of the partial differential equation. It furthermore in- 
dicates that the nonuniqueness problem will occur even 
for very weak shock waves if the airfoil thickness r is 
made sufficiently small. The fact that the Euler results 
do not collapse into a single curve in the similarity plane 
is not unexpected, since the similarity law is only valid 
for the small-disturbance potential equation. The rapid 
drop in dci/da shown in figure 16 for x < 0-8 implies 
a drop in lift, which has been explained as an inviscid 
phenomenon associated with the inefficient production 
of lift that occurs when the upper surface shock wave 
moves to the trailing edge (ref. 17 and ref. 18, p. 656). 

Classically, potential theory is accepted as a good 
approximation to the rotational inviscid flow in the 
transonic range if the shock is weak. The shock strength 
is measured by M % - 1, where M s is the Mach number 
upstream of the shock. The rationale is based on 
the fact that the entropy produced at the shock is 
proportional to the cube of the shock strength. Thus, 
the potential approximation is accepted if 

(M s 2 - l) 3 < 1 (15) 

In practical applications, M s is required to be less 
than 1.3. However, the results of figure 16 are not 
consistent with this requirement. Two potential flows 
are topologically equivalent for affine airfoil profiles 
if they have the same transonic similarity parameter 
X. For example, flow past an NACA 0012 airfoil at 


Moo = 0.820, a = 0°, and flow past an NACA 0006 
airfoil at Moo = 0.879, a = 0° both correspond to 
the same value of the transonic similarity parameter 
X- This value corresponds to conditions on the verge of 
nonuniqueness. It is reasonable to expect that if for one 
of these flow fields the potential assumption ceases to be 
valid, then it should cease to be valid for the other flow 
field as well. The computed surface Mach numbers for 
these two cases are illustrated in figures 17 and 18. It is 
evident from figure 18 that the thinner airfoil does not 
violate the classical criterion. This should be expected, 
since relation (15) is not scaled according to similarity 
rules. The problem can be overcome if we do not allow 
the similarity-scaled entropy 


to be greater than the value it reaches when the 
nonuniqueness occurs. For example, for an NACA 0012 
airfoil near a = 0° , nonuniqueness is observed at a value 
of M s of approximately 1.3; with this information, we 
can compute from relation (16) a value of S which we 
will designate S 0 . For any other affine profile of thick- 
ness r, we must require that the Mach number at the 
shock obey the constraint 

Ms < {1 + S^T 2 / 3 ) 1 / 2 (17) 

For thinner airfoils, this constraint severely limits the 
transonic region for which the potential theory repre- 
sents a valid approximation to an inviscid rotational 
flow. How to extend the constraint given by rela- 
tion (17) or develop some other rule that will prevent 
the application of the conservative potential formula- 
tion in the anomalous range for airfoils of different fam- 
ilies is the subject of an ongoing effort. 

The conservative potential model develops some 
problems when applied to quasi-one-dimensional noz- 
zle flows. It is believed that a heuristic explanation of 
the anomaly can be reached by examining this prob- 
lem. With the nozzle choked, two isentropic solutions 
are possible — one corresponding to subsonic flow, the 
other corresponding to supersonic flow. Of these two, 
the one that appears is determined by matching the 
downstream imposed pressure. If a downstream value 
of pressure is prescribed between these two extremes, 
a shock wave forms in the divergent part of the noz- 
zle or outside the nozzle. The position of the shock 
is established by matching the prescribed downstream 
pressure. The above description is properly modeled by 
the Euler equations. What happens if the problem is 
described by the conservative potential equation? If the 
downstream value of pressure corresponds to the sub- 
sonic isentropic solution, a solution with supersonic flow 
downstream of the throat followed by an isentropic com- 
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pression shock is still possible, since this solution will 
satisfy the downstream boundary condition. Moreover, 
the shock can be placed anywhere downstream of the 
throat and still satisfy the differential equation and the 
boundary condition. In other words, there is an uncer- 
tainty with regard to the shock position. The problem 
comes about because the exact isentropic jump is being 
satisfied and this jump connects the subsonic and super- 
sonic isentropic branches. This is not the case with the 
Euler equations because of the entropy produced at the 
shock, or with the nonconservative potential formula- 
tion because the exact isentropic jump is not satisfied. 
It is believed that this inherent weakness of the con- 
servative formulation to fix the shock position in the 
quasi-one-dimensional problem is also responsible for 
the anomaly observed in two dimensions. 

It can be argued, however, that the downstream 
far-field boundary condition imposed on the two- 
dimensional airfoil problem is specified in terms of 0, 
not (j) x (specifying the back pressure is equivalent to 
specifying (j> x ), and that by specifying <j> the problem 
should then be well posed. However, this argument fails 
to recognize that the shock-wave position is more sen- 
sitive to the pressure level that develops at the trailing 
edge than to the pressure level on the downstream far 
field. Thus, the analogy with the quasi-one-dimensional 
problems remains valid. It is interesting that simi- 
lar ideas for internal flow problems were independently 
developed by Deconinck and Hirsch in a recent paper 
(ref. 19). 

If the nonuniqueness is indeed a problem in the 
shock jump representation, then it could be resolved by 
abandoning the isentropic, mass-conserving shock-wave 
formulation for some other shock representation. (See, 
for example, ref. 20.) However, it seems that to retain 
a rational formulation, some form of shock-fitting will 
be required. 

Concluding Remarks 

The nonuniqueness problem occurring at transonic 
speeds with the conservative potential equation has 
been investigated numerically. Additional evidence 
has been provided which supports the thesis that this 
nonuniqueness is a problem inherent in the conserva- 
tive potential differential equation rather than a con- 
sequence of the finite-difference approximation. The 
hysteresis effect previously observed has been shown to 
be a numerical result of the multivalue nature of the 
computed lift in a certain range of angles of attack at a 
given free-stream Mach number. Since none of the mul- 
tiple solutions obtained with the conservative potential 
formulation seem relevant to the physical problem, a 
breakdown of the theory is indicated. It appears that 
in order to avoid this anomaly, the present conservative 


formulation must be abandoned. A more restrictive cri- 
terion for the validity of potential theory in the tran- 
sonic range has been proposed. This criterion shows 
that shock-wave strength should be measured relative 
to airfoil thickness. The numerical results obtained with 
the Euler code indicate that this nonuniqueness prob- 
lem is not inherent in the inviscid solution but results 
from the approximate treatment of shock waves inher- 
ent in the conservative potential model. However, it is 
felt that more research is necessary to settle this issue 
conclusively. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
October 18, 1984 
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a, deg 


Figure 1. Lift curve obtained with FLO 36 for 
NACA 0012 airfoil with hysteresis effect as in ref- 
erence 3. Mqo = 0.83. 



a, deg 

Figure 2. Lift curve obtained with FLO 36 for NACA 
0012 airfoil at = 0.83. 



Figure 3. Lift curve obtained with FLO 36 for RAE 
2822 airfoil at = 0.75. 
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Figure 7. Convergence histories. NACA 0012 airfoil; M 0 0 = 0.30; a = 10°; 120 x 34 mesh. 




(a) Computed with FLO 36. c; = 1.2733; c<* = 0.0001; 
c m = —0.0100; 192 X 32 mesh; 25 cycles, residual of 
0.452 X 10~ 12 . 
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(b) Computed with FLO 52-S. c t = 1.2613; c d = 0.0026; 
c m = -0.0104; 120 X 34 mesh, 2000 cycles. 


Figure 8. Pressure coefficient for NACA 0012 airfoil at = 0.30 and a = 10°. 
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(a) Computed with FLO 36. c, = 0.0000; c d = 0.0186; 
c m — 0.0000; 192 X 32 mesh; 71 cycles; residual of 
0.660 x 10 -12 . 


-1. 6 1 — 
- 1. 2 - 





.8 



(b) Computed with FLO 52-S. q = 0.0000; c d = 0.0190; 
c m = 0.0000; 120 X 34 mesh; 2000 cycles; residual of ' 
0.165 x KT 6 . 


Figure 11. Pressure coefficient for 


NACA 0012 airfoil at = 0.82 and a = 0°. 




(a) Computed with FLO 36. q = 0.0000; Cd = 0.0854; 
c m = 0.0000; 192 x 32 mesh; 98 cycles; residual of 
0.822 x 10" 12 


(b) Computed with FLO 52-S. c/ — 0.0000; Cd = 0.0577; 
c m — 0.0000; 120 x 34 mesh; 2000 cycles; residual of 
0.124 x 10“ 6 . 


Figure 12. Pressure coefficient for NACA 0012 airfoil at M 0 0 = 0.86 and a = 0°. 
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Figure 13. Shock jump conditions computed with FLO 
36, FLO 52-S, and FLO 12 compared with exact 
jump. 
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Figure 14. Lift curve obtained with FLO 36, FLO 52-S, 
and FLO 12 for NACA 0012 airfoil at = 0.83. 
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Figure 16. Angle made by lift curve in terms of tran- 
sonic similarity parameters for two affine profiles for 
conservative potential and Euler formulations. 
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Figure 17. Computed surface Mach number for NACA 
0012 airfoil at M a 0 — 0.820 and a — 0° with FLO 36. 
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Figure 18. Computed surface Mach number for NACA 
0006 airfoil at — 0.879 and a = 0°, with 
FLO 36. 
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